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Abstract 

A method for finding minimum energy paths of transitions in magnetic systems is presented and 
used to determine the mechanism and estimate the activation energy of skyrmion and antivortex 
annihiliation in nano-systems. The path is optimized with respect to orientation of the magnetic 
vectors while their magnitudes are fixed or obtained from separate calculations. The curvature 
of the configuration space is taken into account by; (1) using geodesics to evaluate distances 
and displacements of the system during the optimization, and (2) projecting the path tangent 
and the magnetic force on the tangent space of the manifold dehned by all possible orientations 
of the magnetic vectors. The method, named geodesic nudged elastic band (GNEB), and its 
implementation are illustrated with calculations of complex transitions involving annihilation and 
creation of skyrmion and antivortex states. The lifetime of the latter was determined within 
harmonic transition state theory using a noncollinear extension of the Alexander-Anderson model. 
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I. INTRODUCTION 


The assessment of the stability of magnetic states with respect to thermal fluctuations is 
an important problem in the theory of magnetism. The preparation of a magnetic system 
in a particular state can be destroyed by thermally activated transitions to other available 
states. Thermal activation also needs to be taken into account when assessing the stability 
of a system with respect to external perturbations such as a magnetic held, contributing, 
for example, to the temperature dependence of hysteresis loops p. 

Thermal stability is a particularly important issue in the context of novel information 
storage devices. Magnetic skyrmions, for example, are localized, noncollinear spin conhg- 
urations demonstrating particle-like behavior and are promising candidates for future data 
storage due to their small size and topological stability. They have been identihed both in 
bulk magnetic materials [2] and in thin magnetic dims [3] . The stability of skyrmions against 
thermal huctuations is an essential prerequisite for their use in memory devices. Originally, 
skyrmions were introduced in the context of elementary particles as conhgurations of contin¬ 
uous fields with topological charge. In a continuum limit, topological protection makes them 
stable against arbitrarily large huctuations. However, in real systems with magnetic mo¬ 
ments localized on atoms, the topological protection is not strict. Instead, skyrmion states 
are separated from the topologically simple (e.g. ferromagnetic) states by hnite activation 
barriers which dehne their stability and lifetime. Vortices and antivortices are other non¬ 
collinear spin conhgurations with potential applications in magnetic information technology 
and spintronics. A methodology for accurately determining the mechanism and rate of an¬ 
nihilation of such noncollinear magnetic states is needed, for example, in the exploitation of 
these kinds of states in devices. 

Thermally activated magnetic transitions are typically rare events on the time scale of 
oscillations of the magnetic moments, making direct simulations of spin dynamics an im¬ 
practical way to calculate transition rates. This separation of time scales, however, makes it 
possible to apply statistical approaches such as transition state theory (TST) [1] or Kramers 
theory |5]. Within the harmonic approximation to TST (HTST) [6] and within Kramers 
theory, the activation energy of a transition, the primary quantity determining thermal sta¬ 
bility, is given by the energy diherence between the local minimum of the energy surface 
corresponding to the initial state and the highest hrst order saddle point (SP) located on 
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a path connecting the initial and hnal state minima. In adaptions of these rate theories 
to magnetic systems [D EHID], the magnitude of the magnetic vectors is either assumed 
to be constant as orientation changes, or it is treated as a fast variable obtained from self- 
consistency calculations for hxed values of the slow variables that specify the orientation m- 
The energy surface of a system of N magnetic moments is then a function of 2N degrees 
of freedom dehning the orientation of the magnetic moments. The description of thermal 
stability of the magnetic states essentially becomes a problem of identifying the SPs on the 
energy surface. For a magnetic system with several degrees of freedom, locating the relevant 
SPs is the most challenging part of the calculation. The difficulty arises from the need to 
minimize the energy with respect to all but one degree of freedom for which a maximization 
should be carried out. The problem is that it is not known a priori which degree of freedom 
should be treated differently. The special degree of freedom for which a maximization should 
be carried out is often referred to as the reaction coordinate. 

An approach which is commonly used in studies of atomic rearrangements but is, however, 
not well justihed, relies on some reaction coordinate chosen before the calculation is carried 
out. For a sequence of values of this coordinate, the energy is minimized with respect to all 
orthogonal degrees of freedom. This method goes by many names and is often reinvented. 
In calculations of atomic systems it has, for example, been referred to as the ‘drag method’ 
because the system is dragged stepwise along the assumed reaction coordinate while it relaxes 
in other degrees of freedom. The maximum energy obtained in this way is then taken to 
be the SP energy. This approach has been applied to studies of magnetization reversal in 
magnetic systems where, for example, a component of the average magnetization vector of 
the system has been chosen as a reaction coordinate [12]. The method works well if the choice 
of reaction coordinate happens to be a good one. It can, however, give inaccurate results if 
the assumed reaction coordinate turns out to be a poor choice. The calculated path traced 
out by this procedure may end up being discontinuous and the maximum energy obtained is 
then not an accurate estimate of the SP energy [13]. In some cases, it is even hard to come 
up with a guess of a reaction coordinate. One example of such a transition in a magnetic 
system is the reversal of a vortex core [TT] which involves formation of a vortex-antivortex 
pair, a complex mechanism for which simple choices of a reaction coordinate would likely 
fail. 

A dehnitive identihcation of the relevant SP for a transition between given initial and 
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final states involves finding the minimum energy path (MEP) between the corresponding 
local minima on the energy surface. Following an MEP means advancing each degree of 
freedom of the system in such a way that the energy is minimal with respect to all degrees 
of freedom perpendicular to the path. An MEP is a path of highest statistical weight 
connecting the initial and hnal states and, thereby, represents a particular mechanism of 
the transition. There can be more than one MEP between given initial and final states, 
each path corresponding to a specihc mechanism and transition rate that can be estimated 
using the HTST or Kramers approaches. The tangent to the MEP provides a good choice of 
a reaction coordinate and the displacement along the MEP represents the advancement of 
the system during the transition from the initial state to the hnal state following a specihc 
mechanism. Any maximum on an MEP is a SP on the energy surface but it is the highest 
maximum that gives an estimate of the activation energy for the corresponding transition 
mechanism. 

The hrst approach we are aware of for hnding MEPs of magnetic transitions involves 
direct minimization of the thermodynamic action of the system [15]. The minimization of 
the action, however, can converge to any path that lies along gradient lines of the energy 
surface and connects the two minima. As a result, an MEPs may not be found but rather a 
path passing through second- or higher-order SPs or even a local maximimum on the energy 
surface. Such ‘false’ paths do not provide an estimate of the thermal activation energy. It 
is expected that the number of such irrelevant paths increases as the number of degrees 
of freedom in the system increases, making it difficult to apply the method to large and 
complex systems. 

The nudged elastic band (NEB) method [131 [16] variations thereof [171120] . are com¬ 
monly used to find MEPs of transitions involving atomic rearrangements. The method 
involves specifying a discrete representation of some initial path between the two minima 
and then using an iterative algorithm to bring the discretization points to the nearest MEP. 
A guess of the reaction coordinate is not needed. Rather, it is determined from the iterative 
optimization of the path. The discretization points represent values of the full set of variables 
of the system and are referred to as ‘images’ of the system. The images provide a discrete 
approximation to a continuous path which should be displaced during the optimization only 
in directions perpendicular to the path [13]. In order to enforce that, a tangent to the path 
needs to be estimated at each image and the magnetic force component in the direction 
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of the tangent removed. In a practical calculation, the number of images is hnite and it is 
important to control the distribution of the images along the path to get adequate resolution 
while keeping the number of images to a minimum and, thereby, minimizing the computa¬ 
tional effort. Good resolution is especially important near the SPs and in regions where the 
path has large curvature. The distribution of images along the path can be controlled by 
adding spring forces between adjacent images. It is important to include only the component 
of the spring force pointing in the direction of the path tangent so as not to interfere with 
the displacement of the path towards the MEP during the optimization. This projection of 
the true force and the spring force using the tangent to the path is referred to as ‘nudging’. 
Without control over the distribution of images, the density of images will drop in regions 
near saddle points where good resolution is particularly important. An alternative way of 
distributing the images is to estimate the length of the path and redistribute the images 
at each iteration using an interpolation between current locations of the images [18]. The 
NEB method has turned out to be a powerful tool for determining transition mechanisms of 
atomic rearrangements, such as chemical reactions and diffusion events. Complex and coun¬ 
terintuitive mechanisms have, for example, been found for catalytic reactions and diffusion 
events controlling the morphology of growing crystal surfaces [221123]. 

The NEB method and variations of it have also been applied to magnetic systems both on 
the atomic scale and on a mesoscopic scale described by micromagnetics [HHIIlEIlElHailSQ]- 
In many cases, these have been rather simple systems with the magnetic vectors rotating 
mainly in a plane (with some exceptions, for example [25l [28l [33|). Several choices of coor¬ 
dinate systems are possible in such calculations. In some cases, spherical polar coordinates 
were used and the length of the magnetic vectors kept constant during the transitions m 
The distance between images in the NEB path was calculated by treating the spherical 
coordinates as if they are linear. This is a rough approximation which leads to reduced 
control of the distribution of the images along the path, especially when some images are 
in the vicinity of the poles. The evaluation of the spring force can also be problematic and 
the springs have in some cases been dropped [21]. But, this is not a recommended practice 
and serves little purpose if the curvature of the conhguration space is properly taken into 
account since the evaluation of the spring force adds only an insignihcant computational 
effort. In other applications, Cartesian coordinates have been used to describe the magnetic 
vectors [31]. This makes the evaluation of the distance between images easier, but the con- 
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tstraint on the magnitude of the magnetic vectors becomes more complicated. It can be 
enforced by projecting the NEB force on the plane perpendicular to the magnetic moments. 
It has been noted, however, that this approach can suffer from convergence problems where, 
for example, the choice of the value of the spring constant becomes critical [3ll |32], unlike 
NEB calculations of atomic rearrangements where the spring constant can be chosen over 
a range spanning many orders of magnitude m ITH] . We show here that this apparent 
sensitivity to the value of the spring constant in calculations of magnetic transitions is not 
present if the tangent to the path is evaluated properly, as described below. 

In this article, an extension of the NEB method to magnetic systems, the geodesic NEB 
(GNEB), is presented and applied to complex transitions involving 3-dimensional rotations of 
the magnetic vectors. The modihcations of the NEB needed for magnetic systems arise from 
the fact that the conhguration space of a magnetic system is a curved manifold due to the 
constraints on the length of the magnetic vectors, which is either assumed to be constant 
or obtained from a separate, self-consistency calculation. As a result, several aspects of 
the NEB method have been modihed, such as (1) the use of geodesics to evaluate the 
distance between adjacent images and to estimate displacements of the images during the 
optimization, and (2) projection of the path tangent as well as the force on the tangent 
space of the manifold of magnetic conhgurations dehned by all possible orientations of the 
magnetic vectors. While various coordinate systems can be used with the GNEB method, we 
have chosen here to use spherical polar coordinates. A special treatment of the coordinates 
is then needed close to the poles. Also, while various minimization methods can be used in 
conjunction with the GNEB method, we present here an algorithm based on an equation of 
motion, analogous to a robust minimization algorithm that has been used extensively with 
the NEB method [13]. The GNEB method can be applied to complex magnetic transitions, 
as is demonstrated here in calculations of skyrmion and antivortex annihilation. The value 
chosen for the spring constant turns out not to be critical, a wide range of values can be 
used, analogous to the NEB method for atomic rearrangements. 

The article is organized as follows. A review of the NEB method for atomic transitions 
is given in Sec. II. The GNEB method is then presented in Sec. III. Results of applications 
of the GNEB method are given in Sec. IV, hrst using a single spin test system to illustrate 
the importance of properly projecting the path tangent as well as problems that can arise if 
springs are not included. Then, the application to skyrmion and antivortex annihilation is 
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described. A discussion and summary is presented in Sec. V. Seven Appendices give addi¬ 
tional technical information: (A) equations for evaluating the path tangent; (B) equations 
for evaluating the geodesic length; (C) detailed description of the Hamiltonian used in the 
single spin test problems; (D) an optimization algorithm based on velocity projection; (E) a 
method for interpolating the energy using a discrete representation of a path; (F) a method 
for verifying that a hrst order saddle point has been identihed; and (G) a brief, step-by-step 
summary of a GNEB calculation. 


II. NEB METHOD 

The basic equations of the NEB are reviewed here in order to make the subsequent section 
on the extensions made in the GNEB method clearer. The NEB method was developed to 
hnd MEPs of transitions where coordinates of atoms change in, for example, chemical reac¬ 
tions, diffusion events or changes in molecular conhguration. The location of each particle 
is given by a 3-dimensional vector, A, and a system of N interacting particles is desribed 
by a 3A^-dimensional vector in Euclidean space R = (fq, r 2 ,... rjv). The fact that the con- 
hgurational space is Euclidean is used explicitly in the calculation of the distance between 
images, to estimate the local tangent to the path and to project the forces. 

Gonsider a chain of Q images, R^, ..., R^~\ , where the endpoints are hxed and given 
by the initial and hnal conhgurations of the transition, but Q — 2 intermediate images 
R’" = ... ,f^), where u = 2,... ,Q — 1, give a discrete representation of a path. 

The position of the intermediate images needs to be adjusted in order to converge on the 
MEP. This is accomplished by displacing the images along the force acting on them so as 
to zero the force. The key idea of the NEB method is the force projection - the nudging. 
It is based on the observation that only transverse displacements of the path should be 
included. A displacement of the images along the path is simply a discrete analog of a 
reparametrization of a continuous path [T3] and does not affect the location of the path in 
conhguration space. Any convenient distribution of the images along the path can be used as 
long as the resolution is high enough. As for any numerical method involving discretization, 
the number of images needs to be large enough to reach convergence. Some reagions along 
the path may require higher resolution than others. One can, for example, choose to increase 
the density of images in regions where the energy is high or the curvature is large HT]. 
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The true force is the negative of the energy gradient, F = —'WE. The transverse dis¬ 
placements are generated by the component of the true force perpendicular to the path. 
Therefore, only this component should be included in the force when the images are dis¬ 
placed. Otherwise, the true force would move the images along the path towards regions of 
low energy and the information about region(s) near saddle point(s), the most important 
part(s) of the path, would be lost. Given an estimate of the unit tangent to the path at each 
image, r'", the perpendicular component of the energy gradient is obtained by subtracting 
the parallel component: 

WE {R') U = WE [R] - {WE {R) ■ f^) rG (1) 

When the distribution of the images along the path is controlled by including springs 
between adjacent images, the spring force has to be projected on the tangent to the path so 
as not to affect the transverse displacements of the images. The parallel component of the 
spring force can be evaluated as 

F;||| = k[\R+^ - R\ - \R - R-^\]f'', (2) 

where k is a spring constant and \R~^^ — R\ denotes the Euclidean distance between 
adjacent images u + 1 and u. 

The total NEB force is then evaluated as: 

F" = -VE(i^")U + F;|||. (3) 

In ordert to form the force projections, it is essential to have a good estimate of the 
tangent to the path at each image, t'^. The simplest approach is to use the line segment 
connecting the previous and later image in the path, i >+1 and u — 1. This, however, has been 
found to lead to instabilities in some cases, slowing down or even preventing convergence 0. 
A better way is to use the neighboring image that has higher energy in either a forward or 
backward hnite difference scheme, and linearly switching between the two schemes when the 
energy of image u is not in between the energy of its two neighbors [35] . 

Some initial path is needed to start an NEB calculation. This can be done in many 
different ways, but the most commonly used method is to generate a linear interpolation 
between the endpoints, = Ri + {i' — 1){Rq — Ri)/{Q — 1). A better starting configuration 
can be obtained by carrying out a procedure based on an interpolation of pairwise distances 


between atoms, the image dependent pair potential (IDPP) method [36], since it prevents 
atoms from coming too close together and generally gives a path that is closer to an MEP. 

When two or more MEPs connect the same initial and final conhgurations, the opti¬ 
mization procednre will most likely lead to convergence to the MEP closest to the initial 
path. In order to hnd the optimal MEP in such a situation, some sampling of the various 
MEPs needs to be carried out. Typically, an NEB calculations involves on the order of ten 
images. When complex paths involving intermediate minima are calculated, it is advisable 
to divide the path at an intermediate minimum and carry out separate NEB calculation of 
each segment of the MEP. 

When the energy along a path is inspected, it is important to use information contained 
in the force along the path in the interpolation scheme as well as the energy at each image. 
This makes it easier to identify, for example, intermediate minima that might be present. An 
interpolation procedure that has proved to be useful in calculations of atomic rearrangements 
is presented in reference |35j . 

Typically, the most important result sought from an NEB calculation is the highest energy 
along the MEP since this gives an estimate of the activation energy for the transition. The 
saddle point will, however, most likely lie in between NEB images and the estimate of the 
maximum energy then be subject to errors in the interpolation between images. In order 
to determine the maximum energy accurately, the highest energy image can be treated 
separately during the iterative NEB optimization and made to move uphill in energy along 
the path. The force on this climbing image (Cl) is calculated by zeroing the spring force 
acting on it and inverting the parallel component of the true force m 

F’^ci = + 2(VE ("4^ 

In this way, the climbing image is made to move uphill in energy along the path but downhill 
in energy perpendicular to the path. The adjacent images effectively dehne the reaction 
coordinate and, thereby, the appropriate climbing direction to reach a hrst order saddle 
point. After the CI-NEB calculation has converged, the position of the climbing image 
coincides with the highest SP along the MEP and gives an accurate value of the saddle 
point energy. 
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III. GEODESIC NEB METHOD 


The state of a magnetic system consisting of N magnetic moments is in principle spec¬ 
ified by 3N parameters - the components of the magnetic vectors in 3-dimensional space. 
However, changes in the magnitude of the magnetic momentum associated with each mag¬ 
netic atom are in most materials much faster than changes in orientation |37|. Therefore, 
the local magnetization and, thereby, the energy of the system can usually be described as 
a function of the orientation of the magnetic vectors only. This is analogous to the Born- 
Oppenheimer approximation in atomic systems where the electronic degrees of freedom are 
assumed to be fast compared with the slowly varying positions of the nuclei. The total 
energy of the system can then be expressed as a function of only the slow degrees of free¬ 
dom. As a result, the configurational space of a system of N magnetic moments contains 
N constraints on the magnitude of the vectors. The NEB method can still be applied in 
such cases, but the constraints would be violated in each iteration and need to be restored 
in some way. Such a ’constrained’ NEB method can suffer from slow convergence similar to 
what has been reported in analogous constrained optimization problems (see, for example, 
ref. [38] and references therein). An alternative is to formulate the calculation as an uncon¬ 
strained optimizattion in properly chosen configuration space. For magnetic systems, instead 
of 3A^-dimensional Euclidean space with N constraints, one can choose a 2A^-dimensional 
Riemannian manifold, TZ, corresponding to the direct product of N 2-dimensional spheres: 

N 

K = ns?. (5) 

i=l 

where Sf is a 2-dimensional unit sphere associated with the Ath magnetic momentum vector. 
The basic equations used in the NEB method need to be revised for this case because the Rie¬ 
mannian manifold is a curved space. The GNEB method presented here is an unconstrained 
optimization method within the 7?.-manifold. 

Similar to the NEB method, the GNEB method involves a chain of images of the system 
giving a discrete representation of a path between the local energy minima corresponding to 
the initial and final states. Adjacent images are connected with springs in order to ensure 
continuity of the path, as in the NEB method. At each image, a local tangent to the path 
needs to be estimated and the force guiding the images towards the nearest MEP is defined 
as the sum of the transverse component of negative energy gradient, i.e. the true force, plus 
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the component of the spring force along the tangent. The position of intermediate images 
is then adjusted so as to zero the GNEB force. 

A particularly important, but non-trivial aspect of the GNEB method is the projection 
of the tangent on the tangent space of the Riemannian manifold. Without it, the true force 
and the spring force interfere with each other leading to uncontrolled behavior of the path, 
as explained below. The tangent space to the 7?.-manifold can be illustrated as follows. 
Gonsider a single spin system where the state is described by a unit vector rh corresponding 
to a point on a 2D unit sphere, S"^. The tangent space to at a point m E is 

dehned as the span of all vectors X perpendicular to rh: 

TrnS^ = {X: X ■rh = 0}. (6) 

For a system of N spins, the tangent space 7^7?. to the 7?.-manifold at a point m = 
(mi, m 2 ,..., mjv) is a direct product of tangent spaces to the unit spheres Sf associated 
with each unit vector m* specifying the orientation of the Tth magnetic moment: 

N 

= (7) 

i=l 

A projection V-j-A of any SWdimensional vector A = (Ai, A 2 ,..., Aat) on the tangent 
space TmT^ can be obtained by subtracting out the component parallel to m* from each 
3 D-vector Ap 

Ai,r = A - {A ■ rhi)rhi. (8) 

VrA is then a composition of Ai^r,..., 

VrA = {A,t, A,t, ■ ■ ■, ^n,t) ■ ( 9 ) 

With the dehnitions given in Eqs. the algorithm involved in the GNEB method 

can now be described. A chain of Q images is constructed, ..., , where 

the endpoints are hxed and given by the local minima corresponding to the initial and 
hnal states, while the location of the Q — 2 intermediate images M'^ = (mj',m 2 , • • • 

V = 2,... ,(5 — 1, is adjusted so as to zero the GNEB force acting on them. The basic 

force is dehned in a way analogous to the NEB force, as in Eqn.(|^, with E{Ry) replaced 
by But, now the tangent used to project the true force and the spring force has to 

be projected on the tangent space of the 7^-manifold. First, the tangent is estimated in a 
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way analogous to the NEB method, see Appendix]^ and then it is projected on the tangent 
space 


ff — 


Vtt^ 


( 10 ) 


The perpendicular component of the energy gradient can then be obtained by subtracting 
out the parallel component 


S/E {M^) U = S/E {M^) - (VE (M 


T-y) T-y. 


( 11 ) 


The parallel component of the spring force is evaluated using geodesic distances 


F;||| =k{L - L T- 


V 

r* 


( 12 ) 


Here, L (M''+', M") and L (M", M"-') are geodesic distances between images z/ + l, z/ and 
z/, 1 / — 1, respectively. Since the 7?.-manifold is equal to a direct product of (see Eq. 
the distance between points v and p is given by: 


L (M", Af") = . 


(13) 


Here are great-circle distances between points v and /r on the Tth unit sphere. can be 
computed using, for example, spherical law of cosines, haversine or Vincenty’s formulae [SU^ 
53] (see Appendix 0. Here, the same spring constant, zt, is used in the applications of 
the method, but a non-uniform distribution of the images along the path can be obtained 
by choosing different values of the spring constant for each geodesic segment, using stiffer 
springs in regions where higher resolution of the path is desired. 

The force F'^ obtained from Eq. (|^ by inserting VE (M’^) |x from Eq. (pTj) and from 


Eq. (12) is in general not in the tangent space of TZ even though the tangent, G Tm''TZ, 
is. The force that should be applied to the images in the GNEB method, therefore needs to 
be projected on the tangent space to satisfy the magnetic constraints: 


fS«EB = Pr(-Vi5(M-)U + F;|||) (14) 

where the projection, Vr, is calculated using Eqs. ([^ and (|^. The GNEB optimization 
procedure iteratively zeros the force to bring the images to the nearest MEP. 

The key difference between the tangent used in the GNEB method, Tf, and the tangent 
used in the NEB method, r^, is that the former lies in the tangent space of the 7?.-manifold 
of image z/, see Fig. la. This is needed for the true force to be decoupled from the spring 
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force. If the NEB tangent is nsed for magnetic transitions, withont projecting on the tangent 
space, interference occnrs between the trne force and the spring force and, thns, nncontrolled 
behavior of the images. The reason for this is illnstrated in Fig. [^for a single spin system. 
Three neighboring images, z/ — 1, n and z/ + 1 are placed eqnidistantly along a geodesic 
path which lies in the plane of the hgnre. The forward-difference approximation is nsed to 
dehne the local tangent to the path at image z/: — fhy) / — m^|. Since 

the distance between adjacent images is the same in this case, the spring force is zero. 
Even if the transverse component of the trne force is calcnlated by snbtracting ont the 
component parallel to the NEB tangent and projecting on the tangent space so as to 
satisfy the magnetic constraint, the resnlting force still has a component along the path. 
The trne force then affects the distribntion of images along the path. This will canse image 
V to move towards the lower energy region of the path, i.e. slide down. Control of the 
arrangement images along the path is then lost and, as a resnlt, the resolntion of the path 
in the critical regions becomes poor. The valne chosen for the spring constant then also 
becomes an important issne. This behavior is demonstrated in an application to a single 
spin test problem in the following section. The problem can be avoided by projecting the 
tangent on the tangent space as illnstrated in Fig. la. 

An initial path is needed to start the GNEB calcnlation. In many cases it is snfficient to 
place the images initially along the geodesic path between and . This is analogons 
to the straight-line interpolation between the initial and hnal point often nsed to start NEB 
calcnlations of atomic rearrangements. A simple way to obtain the intermediate images 
= ■ ■ ■ 1 ^'n) along the geodesic path is to rotate each magnetic momentnm 

from the initial orientation to the final one using Rodrigues’ rotation formula: 

mj' = ml cosozf -1- {ki x m[) sincnj". (15) 


Here m\ is the i-th unit magnetic vector in image z/th and uj\ = {u — IjAcui is an angle of 
rotation. If Q images are used in the chain, then Aui = uJi/{Q — 1), where Ui is an angle 
between vectors mf and rhf. A unit vector ki describing an axis of rotation is dehned as: 


(mf X mf) 
{ml X mf) 


(16) 


As with the NEB for atomic rearrangements, when multiple MEPs are present, the opti¬ 
mization procedure will most likely lead to convergence to the MEP closest to the initial 
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FIG. 1. Illustration of the force projection used in the GNEB method and the error arising if 
the tangent is not projected on the tangent space. Three equidistant images, u — v and + 1, 
marked with filled circles, are located along a geodesic path lying in the plane of the figure. The 
NEB tangent to the path at image u, tneb, is defined as a unit vector pointing to a higher energy 
image v + 1. (a) The tangent projected on the tangent space, tgneb = "tT) is used to project out 
the parallel component of the gradient, VE, to obtain the transverse component, VEj^, which does 
not move the image along the geodesic path. The distribution of images along the path is, thereby, 
not affected, (b) A problem occurs if the tangent is not projected on the tangent space and tneb 
is used to estimate the transverse component of VE. A non-zero component of the force acting on 
the image remains in the tangent space, as shown with a red arrow. The true force then interferes 
with the distribution of images along the path. 

path and some sampling of the various MEPs is needed to hnd the optimal one, or - more 
generally - all the relevant ones. 

A summary of the GNEB algorithm is given in AppendixA climbing image can be used 
in a GNEB calculation, a GI-GNEB method, analogous to the GI-NEB method described in 
the previous section. The highest energy image is chosen to be the climbing image and the 
force on it is only the inverted parallel component of the true force without the addition of 
a spring force 

= Vt{-VE + 2{VE ( 17 ) 

The GI-GNEB method is most useful when the energy barrier is narrow and sharp compared 
with the length of the path. 

As for NEB calculations of atomic rearrangements, it is important to interpolate the 
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energy along GNEB paths using information contained in the force along the path as well 
as the energy at each image, as described in Appendix This is important for obtaining 
an accurate estimate of the saddle point energy as well as to identify possible intermediate 
minima. 

The formalism described above for the (CI-)GNEB method is independent of the coor¬ 
dinate system. We have chosen here to use spherical polar coordinates since they naturally 
fulhll the magnetic constraints. A conhguration of a system of N spins is specihed by the 
coordinates ( 6 *i, 0 i, 6 * 2 , 02 , • • • 6 *^, 0 Af) and all displacements are then automatically within 
the 7^-manifold. 


IV. RESULTS 

A. Test problems 

The GNEB method is hrst illustrated with calculations of the MEP of a transition in¬ 
volving a single spin where the energy surface can be visualized easily. The energy is given 
by a sum of Gaussian functions of geodesic distances, as described in Appendix 

First, several minimization calculations were carried out where parameters of the energy 
surface and starting points were chosen randomly. The minimization algorithm is based on 
an equation of motion where the velocity is projected in such a way as to converge on a 
minimum, as described in Appendix]^ The test calculations were done to assess the robust¬ 
ness of the minimization algorithm, especially near the poles. In some of the calculations, 
the energy surface was rotated so that one of the minima conicided with a pole. With the 
special treatment of the polar regions, described in Appendix the minimization converged 
witout problems. Furthermore, the number of iterations needed to reach convergence did 
not tend to be larger when the minimum coincided with a pole. 

Fig. g shows a GNEB calculation of the MEP for a transition of the spin from one 
minimum to another. The initial distribution of images was generated along the short 
geodesic path connecting the minima. The geodesic path is obviously not an optimal one 
since it lies close to an energy maximum. The images where then iteratively brought to the 
MEP using the GNEB method with the minimization algorithm described in Appendix 
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FIG. 2. A GNEB calculation of an MEP in a single spin system. The positions of the images are 
shown with filled circles. The initial path was chosen to lie along the short geodesic path between 
the energy minima (indicated by blue color) and turned ont to lie near a maximum (indicated by 
red color). The converged path is shown in pink. Two intermediate configurations of the path are 
also shown. The saddle point is marked with a cross. Q = 15 images were used in the calculation. 

1. Effect of not projecting the tangent on the tangent space of TZ 

As illustrated in Fig. [T| the true force and the spring force are not decoupled properly 
unless the tangent to the path is projected on the tangent space. Without this projection, 
interference between the forces lead to loss of control of the images, such as sliding down 
from the saddle point region towards the minima. Fig. |^a) shows the results of such a 
calculation. After 200 iterations the images have moved to the vicinity of the MEP, but 
the spacing between images is uneven and the resolution of the path in the barrier region is 
unsatisfactory. Since the parallel component of the true force is not projected out completely, 
the images slowly slide down from the barrier region. After 400 steps most of the images 
have slid down into the potential wells and the information about the saddle point region is 
lost. It is still possible to obtain a reasonable estimate of the MEP without the projection 
on the tangent space but it requires choosing just the right value for the spring constant. 
Such behavior has indeed been reported before [32]. 

When the tangent to the path is correctly defined by including the projection on the 
tangent space, the method is, however, stable and quite insensitive to the choice of the value 
of the spring constant. In fact, the GNEB calculations for the test problem in Fig. [^were 
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found to give nearly identical results as the value of the spring constant was varied over 
hve orders of magnitude. For complex transitions involving multiple spins and rotation in 
3-dimensions, proper convergence to the MEP within a tight force tolerance (such as the 
10“® meV/radian tolerance specihed in the present calculations) may not occur for any value 
of the spring constant unless the tangent to the path is projected on the tangent space. 




FIG. 3. Same test problem as in Fig. The images are shown with filled circles and the saddle 
point with a cross, (a) Results of calculations where the tangent is not projected on the tangent 
space. The initial distribution of the images was chosen to be along the geodesic path connecting 
the energy minima. The path after 200 time steps is shown with a dashed line and the path after 
400 steps is shown with a solid line. The parallel component of the true force is not completely 
projected out in this case and that makes the images slide down from the barrier region, (b) Results 
of optimization without springs but including the full GNEB force projections. Because of kinks 
on the path, a wiggling motion during the iterative minimization enables the images to gradually 
slide down in energy. As a result, the resolution near the saddle point becomes poor. Small random 
displacements were added to the initial arrangement of the images along the geodesic. 


2. Effect of skipping the springs 

Since the forces are decoupled in the GNEB method, the value of the spring constant 
is not critical. For some simple systems, reasonable results can even be obtained with the 
NEB without the inclusion of springs 125, but control of the distribution of the images is 
then lost and the minimization does not converge properly. Fig. [^b) shows the results for 
the same test problem as the one shown in Fig. The images in the initial path were 


17 









arranged along the geodesic between the two minima with small random displacements 
added before starting a GNEB calcnlation withont springs. The calculation converges in 
this case, i.e. perpendicular component of the true force drops below the set tolerance, but 
the hnal distribution of the images is uneven, as most of the images have slid down towards 
the energy minima. Without the springs, the method becomes unstable with respect to 
small perturbations. The inclusion of the springs as in Fig. ensures convergence to the 
MEP with an equal distribution of images when the same value of the spring constant is 
chosen for all the springs. 


B. Annihilation of a skyrmion 

The GNEB method was used to calculate the MEP for annihilation of a skyrmion to 
form a homogeneous, collinear state. While these are topologically distinct magnetic states, 
there is a hnite energy barrier for transitions between the two. The MEP not only gives an 
estimate of the activation energy barrier but also reveals the microscopic mechanism of the 
skyrmion annihilation, or - in reverse - skyrmion formation. 

The emergence of magnetic skyrmions is the result of an interplay of multiple interactions, 
where the Dzyaloshinskii-Moriya (DM) interaction and noncollinear isotropic exchange are 
believed to play a crucial role [HI [32] • A discussion of the microscopic origin of skyrmion 
states is, however, outside the scope of this article. Here, skyrmion conhgurations are 
obtained using a simplihed phenomenological model based on a Heisenberg-type Hamiltonian 
including isotropic exchange, DM interaction as well as interaction with an external held. 
The energy of the system is dehned as: 

^ TV ^ N N 

E = --J^rhi-rfij - - ^Dij ■ [miX rhj] ~ B^Mi (18) 

{i-j) (*J> * 

Here, J is the isotropic exchange parameter, the Dij are the DM vectors and B denotes 
an external magnetic held. Values of the parameters of the Hamiltonian were taken from 
Ref. [33] . The exchange interaction is included only between nearest neighbors (indicated 
by angular brackets). The spins are located in the XY-plane at the nodes of a 2-dimensional 
square lattice. The size of the simulated domain was 21 x 21 spins, which is signihcantly 
larger than the size of the skyrmion. Rigid boundary conditions were applied, essentially 
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mimmicking an infinite 2-diniensional layer. The magnetic field was applied along the normal 
to the plane. The DM vectors lie in the plane as shown in Fig. 



FIG. 4. Minimum energy path showing the mechanism and activation energy for the annihilation 
of an isolated skyrmion. The filled circles indicate images along the converged path obtained using 
the CI-GNEB method. The interpolation method described in Appendix |D] was used. Insets show 
the orientation of magnetic moments at three configurations. The color indicates the value of the 
z-component of the magnetic vectors (red -H- up, green -H- down). Only a part of the simulated 
system is shown. The ground state corresponds to collinear ordering of the magnetic vectors. 
The thin, blue arrows in the left-most inset show the direction of the pseudovectors for the DM 
interaction of a central moment with its four nearest neighbors. 

Minimization calculations starting from random initial configurations of the spins revealed 
two stable states in the system. The global energy minimum represents a collinear state 
with all the magnetic vectors aligned along the field B. A noncollinear, metastable state 
corresponding to an isolated skyrmion was also found. The skyrmion center coincides with 
a lattice site, see Fig. The CI-GNEB calculation was started by generating 19 images 
along the geodesic path between these minimum energy configurations and then applying 
the minimization algorithm described in Appendix until the force on each magnetic 
momentum had dropped below 10“® meV/radian. The calculation converged to the MEP 
without any problem. Fig. shows the result. The MEP reveals a mechanism for the 
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annihilation of the skyrmion. At hrst, the skyrmion moves as a whole withont changing its 
size and shape, shifting its center in between the atomic sites. This shift involves almost 
no change in energy. After that, the spins rotate symmetrically cansing the skyrmion to 
gradnally shrink and eventnally disappear. This is a challenging calcnlation involving ont- 
of-plane rotation of the magnetic momenta, as can be seen from the insets of Fig. 

The activation energy for the annihilation of the skyrmion state was calcnlated to be 40 
meV, bnt 50 meV for its formation from the collinear state. The saddle point conhgnration 
is shown in an inset of Fig. There, the fonr central spins lie in the XY-plane and point 
towards each other. Each of them is perpendicnlar to the nearest neighbors bnt antiparallel 
to next nearest neighbors. 


C. Annihilation of an antivortex 

The noncollinar Alexander-Anderson (NCAA) model has been fonnd to describe well 
magnetism of itinerant electrons in 3d transition metals |l5l 06]. While the evalnation 
of the energy reqnires self-consistency calcnlations, analytical forces provided by a force 
theorem m make it relatively easy to navigate on the energy surface and hnd local minima 
corresponding to (meta)stable magnetic states with possibly complex, non-collinear ordering 
of the magnetic moments. Calculations of a 7x7 atomic row island of Fe atoms on a substrate 
have been shown to lead to an antivortex state when the parameter values are slightly 
different from the optimal values of an Fe crystal under ambient conditions m- Such a slight 
change in the model parameters could be the result of an external perturbation such as an 
external electrical held or the presence of impurities or defects. The metastable antivortex 
state was found by starting from a random initial orientation of the magnetic moments and 
then minimizing the energy to bring the system to a minimum energy conhgnration. This 
state can be described as a symmetrical, saddlelike arrangement of the magnetic moments 
in the center of the island. The total in-plane magnetization is zero, while the out-of-plane 
magnetization is non-zero mainly due to four magnetic moments near the center of the island 
pointing out of plane. 

The CI-GNEB method was used to calculate the MEP for the transition between the 
antivortex state and the homogeneous ground state, and the results are shown in Fig. 
The calculation was considered converged when the force on each magnetic momentum 
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FIG. 5. Minimum energy path showing the mechanism and activation energy for the annihilation 
of an antivortex state obtained in a noncollinear Alexander-Anderson model of a 7x7 island of Fe 
atoms. The filled circles indicate images along the converged path obtained using the CI-GNEB 
method. Insets show the orientation of the magnetic moments at various points along the path, 
as indicated by arrows, colored with the same scheme as in Fig. The ground state corresponds 
to collinear ordering of the magnetic vectors. The interpolation method described in Appendix [D] 
was used. 
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had dropped below 10“® meV/radian. From these results, the activation energy for the 
annihilation as well as the formation of the antivortex state can be obtained. As the system 
progresses along the MEP, the saddle-like excitation moves along the diagonal of the island 
towards one of the corners (the upper right hand corner in the insets of Fig. where it leaves 
the island. As with the skyrmion annihilation, Ending this MEP is a challenging calculation 
involving out-of-plane rotation of the magnetic momentum vectors, as can be seen from 
the insets of Fig. The maximum energy configuration was verified as a SP configuration 


using the method described in Appendix G The activation energy for the annihilation 
of the antivortex state was calculated to be 10 meV, but 55 meV for its formation from 
the collinear state. The pre-exponential factor in the rate constant was determined using 
HTST and found to be 1.4x10^^ s“^ for the annihilation of the antivortex. For the reverse 
transition, the pre-exponential factor was larger, 2.6x10^^ s“^, showing that the vibrational 
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entropy of the antivortex state is larger than that of the homogeneous state. 


V. DISCUSSION 

The results presented here on the annihilation of the skyrmion and antivortex states are 
preliminary in that they have been obtained for only particular systems and Hamiltonians. 
They are presented here mainly to illustrate possible transition mechanisms and give an 
indication of the thermal stability in these simple systems. It is not known at this time 
whether similar transition mechanism will be obtained for other models supporting these 
noncollinear states. This remains to be seen. 

The existence of a local minimum on the energy surface of the NCAA model corresponding 
to an antivortex state demonstrates the complex exchange interactions included in the model 
even though it contains only a few parameters. Such a behavior cannot be obtained within a 
Heisenberg-type model unless additional phenomenological terms and additional parameters 
are introduced in the Hamiltonian. In the NCAA model, complex noncollinear states appear 
quite naturally. It may be possible to obtain skyrmion states within the NCAA model. 
Again, this remains to be seen. 

Most importantly, the calculations of the annihilation of the skyrmion and antivortex 
states presented here illustrate the applicability of the GNEB method to complex magnetic 
transitions. While the NEB method has been used already in several studies of magnetic 
transitions [HHIIl EB ElHSl SO] convergence problems and sensitivity to parameters such 
as the spring constant have been noted and in light of the analysis presented here it is likely 
that rigorous convergence, dehned here as a drop in magentic forces below 10“® meV/radian, 
was in some cases not obtained. The extension of the NEB method presented and used here, 
the GNEB method, resolves these convergence issues and provides a tool for the study of 
more complex magnetic systems than could be studied with the NEB method. 

The GNEB method assumes a certain shape of the conhguration space: a direct product 
of 2D spheres associated with each magnetic moment in the system. This requirement is 
clearly fulhlled by the quasi-classical Heisenberg-type models and micromagnetics, where 
the length of magnetic vectors is usually hxed and does not depend on orientation of the 
spins. It is also fulhlled in the quantum-mechanical treatment of the system within an 
adiabatic approximation, where the magnitude of the magnetic moments is a function of 
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the orientation dH E?]. The magnitnde of the magnetic moments is assumed to adjust 
instantaneously to changes in orientation. The energy surface is entirely dehned by the 
orientation of the magnetic momentum vectors. In general the magnetic force can contain 
a term proportional to the derivative of the energy with respect to the length of magnetic 
moments, which accounts for the change in magnitude of magnetic moments as magnetic 
vectors rotate. However, this term vanishes for all important cases such as density functional 
theory or NCAA model calculations as shown by a magnetic force theorem |IIIS7|. 

Here we present some practical information which we hnd useful for a successful GNEB 
calculation. If interpolation of the energy along the path using the method described in 
Appendix gives an indication of an intermediate minimum, it is recommended to divide 
the path up into segments and study transitions between each pair of adjacent minima 
separately. The location of the minimum should hrst be carefully identihed by carrying out 
an energy minimization and then calculate the MEPs between each pair of adjacent minima 
in separate GNEB calculations. The full MEP is, thereby, obtained as a composition of 
segments of MEPs, each passing through only one SP. 

The question arises when to use GNEB and when to use CI-GNEB. The latter is partic¬ 
ularly efficient when the energy maximum along the path is well dehned. However, there 
are important exampes of magnetic transitions where the barriers are hat, making it hard 
to identify the highest energy image. Flat barriers appear, for example, when the mecha¬ 
nism of a magnetic transition involves domain wall propagation without signihcant change 
in energy In such cases, it is better to use the regular GNEB method rather 

than the GI-GNEB method. A good strategy for hnding MEPs is, therefore, to start the 
calculation using the GNEB method and inspect the energy along the path after a few tens 
of iterations or when the GNEB forces have dropped below some predehned value (signih- 
cantly larger than 10“® meV/radian). If there is an indication of a flat barrier, then continue 
using GNEB until convergence (see Appendix]^. If the path seems to have a well-dehned 
maximum, then switch to the GI-GNEB method. During the GI-GNEB calculations, the 
energy of the climbing image might become lower than that of some other image in the 
chain. In this case, the climbing image should to be reassigned. 

The GNEB method does not imply a particular optimization method. Any optimization 
method based only on the evaluation of the forces, i.e. the hrst derivatives of the energy 
with respect to the conhguration parameters, can be used. The objective function corre- 
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spending to the GNEB forces is not well defined, however, because of the force projections 
and magnetic constraints. It remains to identify the most efficient optimizers for the GNEB 
method, similar to what has been done for the NEB method |1S]. The method presented 
in Appedix D is stable and reliable, but probably not optimal in terms of efficiency. It 
is likely that a switch to some quadratically convergent optimization method is beneheial 
when the images are relatively close to the MEP, while the use of a more stable method, 
such as the one presented in Appedix D, is safer when starting from the arrangement of 
the images along the geodesic path. Since the most important point along the path is the 
highest energy SP, it can be efficient to converge the MEP only to a rather high tolerance 
but then converge with a low tolerance on the SP using, for example, the minimum mode 
following method 150]. 

To summarize, we emphasize that the GNEB method differs from the NEB method in 
that the curvature of the conhguration space of magnetic systems is taken into account. A 
particularly important aspect of the GNEB method is the projection of the tangent to the 
path on the tangent space. If this projection is not implemented, a rigorous convergence of 
to the MEP will likely not be obtained. In fact, we tried to reproduce the results obtained 
with the (GI-)GNEB calculations on the skyrmion and antivortex annihilation using simpler 
implementations of the NEB method for magnetic transitions where projection of the tangent 
to the path on the tangent space is not included (as in refs. mm), but could not reach 
convergence to the MEP for any value of the spring constant. 
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Appendix A: Evaluation of path tangent 


The tangent to the path at image can be estimated using the coordinates of the two 
adjacent images in a straightforward way as 

but this can lead to the formation of kinks in the path and, as a result, cause instabilities 
that slow down or even prevent convergence to the MEP [T3]. This problem can be reduced 
by using a better definition of the tangent based on either forward or backward difference 
depending on the energy of the image [35] : 


{ rjj, if E’^+^ > > E'^-^ 

rjj, if < E'^-^ 


(A2) 


where rjj = and rjj = and E'" = E If image v is close to 

an energy minimum > E'^ < E'^~^) or to an energy maximum < E'^ > E'^~^) 

along the path, the tangent is taken to be a weighted average of rjj and rjj so as to ensure 
a smooth switch between the forward and the backward difference schemes: 


where 


max 


^ IrjjAEj;,, + rfjAE;;,,, if 

[rl^E^rmn + if < E’^-\ 

(l^^+i _ E'^\ , - E^l), and AEj;^^ = min(|E^+^ - E'^\ , \E 


(A3) 


Appendix B: Evaluation of the geodesic distance 

An essential aspect of the GNEB method is to evaluate the geodesic distance between 
two locations of a magnetic vector. Three different equations for this calculation are given 
here. The formulas are mathematically equivalent but they differ in numerical stability. 
The spherical law of cosines [ST] 

= arccos(mj' ■ mf). (Bl) 

is the simplest formula but it can lead to significant round-off errors when the distance 
between the two points on the unit sphere is small. 
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The haversine formula works well for small as well as large distances but is more compli¬ 
cated [52] 


,1/,/i 

'-i 


2 arcsin 




(B2) 

Here, 6^ and (j)^ are polar and azimuthal angles dehning orientation of The only problem 
arises for antipodal points (on directly opposite sides of the sphere), where the formula suffers 
from round-off errors. 

Vincenty’s formula is [53] 


= arctan2(|m^ x mf| ■ mf). 


(B3) 


where the function arctan2(?/, x) is dehned as follows: 


arctan2(?/, x) = < 


arctan —, if x > 0 

X 

y 

arctan —h tt, if ?/ > 0, x < 0 

X 

y 

arctan- tt, if r/ < 0, x < 0 

X 


-f- 
^ 2 ’ 

_ 

2 ’ 

undefined. 


if r/ > 0, X = 0 
if r/ < 0, X = 0 
if 1/ = 0, X = 0. 


(B4) 


It works well for all caes but is the most complicated expression and requires the dehnition 
of the arctan2 function because the usual arctan function gives values in between —7r/2 and 
7r/2 while the geodesic length on the unit sphere should be between 0 and tt. We have used 
Vincenty’s formula in the calculations presented here. 


Appendix C: Summary of GNEB algorithm 

A brief, schematic summary of a GNEB calculation is as follows: 

1. Identify two minima on the energy surface correspondig to the initial and hnal conhg- 
urations, = (m{, • • •, ^n) and = (mf, m^,..., , between which an MEP 

is to be found. 
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2. Create an initial path by constructing a chain of intermediate images between and 
M^. When the geodesic path is used to construct the initial path, an intermediate image, 


I/, with coordinates ... ,m'^) is generated using Eqs.(15) and (16). 

3. For each image, set the velocity, V', to zero. 

This completes the initialization of the iterative algorithm. 

4. Calculate the energy of the images, E {M^). 

5. For each image, calculate the tangent to the path, as described in Appendix]^ Then 


project the tangent on the tangent space, Fqs.(lO) 


6. For each image, calculate the transverse component of the true force, Fq.([II]). 

7. Calculate the geodesic distance between pairs of adjacent images, Fq.( 


8. For each image, calculate the parallel component of the spring force, Fq.(12). 


9. For each image, calculate the total force by adding the transverse component of the 
true force projected on the tangent space to the parallel component of the spring force. 


Fq.(14). 


10. If the GNFB force on any one of the magnetic vectors in any one of the images is larger 
than the tolerance, Ftou then update the position and velocity according to the optimization 
algorithm in Appendix and go back to step 4. Otherwise stop, and interpolate the energy 
along the path using the method in Appendix and verify that the highest energy point is 
a hrst order saddle point using the approach in Appendix 


Appendix D: Interpolation between images 

The hnal, relaxed conhguration of the images obtained from a GNFB calculation gives a 
discrete representation of the MEP. Except for some simple, highly symmetrical cases, the 
highest energy image in the discrete representation of the path will likely not coincide with 
the energy maximum of the MEP, i.e. a SP on the energy surface. In most cases, the SP 
lies in between the images. An accurate interpolation of the path between the images is, 
therefore, important. It is also useful for identifying intermediate minima along the path. 
Not only the energy at each image is available for the interpolation but also the component 
of the energy gradient along the path at each image, = (VE ■ -rf). This input can 

be used to construct a piecewise cubic polynomial interpolation. The energy of the system 
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is approximated as 


E(X) — a,y{\ — + b^{\ — + C;y(A — X^ + d^, X G [Aj^; Aj^+i] (Dl) 

where A is the displacement along the poly-geodesic MEP and X^, is given by 

{ 0, if z/ = 1, 

(D2) 

A,_i + L(M"-i,M"), if z/ = 2,...,g. 

For each segment, the [Ai,;A,y+i], parameters a^, Cy and are dehned so as to enforce 
continuity of the energy and its derivative at both ends of the segment: 


d„ = E' 



2 - E’') 

(D3) 

(Aj/+i -l- Ai,)^ 

(Ai/+i -l- XYf 

^ 26’' 

3 {E’'+^ - E’') 

1 , , o s 

(D4) 

Xu+l + Xy 

(Ajy+l -|- Xy) 

5Y 


(D5) 

E’'. 


(D6) 


The derivatives of the energy, 5^-, provide important information about the path. In some 
cases, intermediate minima on the MEP can be identified which could not be seen when 
only values of the energy are used in the interpolation. When an intermediate minimum is 
found, it is advisable to split the path and calculate the two parts of the MEP separately. 

The orientation of the magnetic vectors at the saddle point can be obtained from the 
value of A at the maximum of the energy curve. 


Appendix E: Energy surface for the test calculations 

For the simple test problem calculations, the energy of the system is taken to be a sum 
of Gaussian functions of the geodesic distances 

^ r 

E(m) = ^ajexp , (El) 

i=i ^ ^ 

where lj{m) is the geodesic distance between the endpoint of vector m and a center of the 
j-th Gaussian, mj. Parameters dehning the center of Gaussians, ruj, their amplitude, aj, 
and width, aj were chosen so as to form multiple minima, maxima and saddle points on the 
energy surface. 
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The values of the parameters of the Gaussian functions, i.e. the amplitude (oj), width 
(uj) and position of the center (in spherical coordinates 6i and (pi) used in calculations shown 
in Figs. 2 and 3 are listed in the table below: 


i 

di 


e^ 

(pi 

1 

-2.500 

0.200 

0.600 

3.700 

2 

-2.000 

0.150 

0.650 

6.300 

3 

4.000 

0.300 

0.150 

5.000 

4 

0.656 

0.080 

1.523 

4.280 

5 

1.609 

0.341 

1.662 

4.899 

6 

0.696 

0.303 

1.753 

3.963 

7 

1.753 

0.488 

1.404 

5.777 

8 

0.429 

0.276 

0.132 

5.655 

9 

1.659 

0.862 

1.831 

3.974 

10 

1.222 

0.568 

1.291 

4.539 

11 

1.302 

0.402 

0.020 

1.968 

12 

1.685 

0.743 

2.494 

1.141 

13 

1.287 

0.800 

1.438 

3.046 

14 

0.249 

0.282 

0.780 

2.869 

15 

0.304 

0.443 

1.609 

0.611 

16 

0.656 

0.430 

1.361 

4.123 

17 

0.057 

0.601 

2.260 

2.854 

18 

1.337 

0.236 

2.245 

2.135 

19 

0.075 

0.813 

0.561 

5.211 

20 

0.652 

0.776 

2.232 

4.565 

21 

0.392 

0.453 

1.741 

4.825 

22 

0.014 

0.565 

0.860 

3.719 

23 

0.111 

0.484 

1.648 

0.490 

24 

0.884 

0.160 

2.867 

5.354 
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Appendix F: Iterative optimization algorithm 


The GNEB force needs to be zeroed by moving down-hill in energy towards an MEP. 
Although many techniques are available for such an optimization problem, standard iterative 
algorithms where the search direction at each step is defined by the force are not efficient for 


curved manifolds such as TZ. Even though the force defined in Eq. (14) is correctly confined 
to the local tangent space, any finite displacement along a straight line in the direction of 
the force brings the system out of the 7^-manifold and the magnetic constraints then need to 
be restored at each iteration. It is preferable to advance the system along a geodesic of the 
curved 7^-manifold [5H1 El] • We have used a simple but efficient method based on an equation 
of motion where the velocity is damped by including only the component in the direction of 
the force. The method is referred to here as velocity projection optimization (VPO). It is 
analogous to a method that has been used in NEB calculations for unconstrained Euclidean 
configuration space [T3|. Here, it is adapted to magnetic systems by advancing the system 
along geodesics, thereby automatically fulfilling the magnetic constraints. The details are 
given in the following subsection. Some special considerations arise in the regions close to 
the poles. A method for resolving those issues is described in the second subsection. 


1. Velocity projection optimization on a sphere 

In the VPO method, the coordinates and velocities are updated according to equations of 
motion. In every iteration, only the component of the velocity parallel to the force is kept, 
unless it is pointing in a direction opposite to the force, at which point it is zeroed. The 
advantage of this formulation is that if the system moves in the right direction, it accelerates 
and approaches faster the region of an optimum. If the algorithm overshoots and the system 
has gone passed the optimum, the velocity is zeroed. 

In order to be able to use the VPO method for magnetic systems, the evolution of the 
system has to satisfy the magnetic constraints. While the Landau-Lifshitz equations describe 
the spin dynamics, they cannot be used in this optimization algorithm since they do not 
involve mass and acceleration. Instead, the motion of a point mass on the manifold TZ is 
used here. This motion is conveniently described in terms of spherical coordinates, 6i and 
(pi, which automatically account for the curvature of the manifold TZ. For a system of N 
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magnetic moments, the following set of 2N coupled equations are used: 


m9i — m(f)l sin 9i cos 9i 

= f" 

J 1 1 

(Fl) 

m(f)i sin 9i + 2m4>f cos 9i 

_ 

Ji ’ 

(F2) 


Here, m is the mass, {Oi,(j)i) are polar and azimuthal angles dehning the orientation of i- 
th magnetic moment, and (/f,//) are projections of the i-th force on the orthogonal unit 
vectors, (ef,ef), in the direction of increasing {9i,(j)i). These equations are similar to the 


equations of motion for a spherical pendulum [55]. Eqs. (FI), (F2) are then transformed to 
coupled hrst order equations for position and velocity: 


Vi = vfvf cot 9i + —/f, 
m 

vf = -v-vf cot 9i + 

I I i 

Oi= y. 



(F3) 

(F4) 

(F5) 

(F6) 


These equations are solved iteratively using a numerical scheme. We have found that even 
the simplest Fuler integrator works well. 

At each time step, the velocity is modihed depending on the sign of its projection on the 
direction of the force. More specihcally, if is a 2Wdimensional vector of the force with 
components /f,..., /^, //,..., /^, and V = (^vf ,..., v%, vf, ..., j is a 2Wdimensional 
velocity vector, then 


(V ■ 5r) 3 ^/ , if (V ■ 5r) > o, 

yrev ^ I ^ 

lo, if(V-9")<0. 

Here, superscript ’rev’ stands for ’revised’. The revised velocity is used in the following 
iteration. If the force is consistently pointing in a similar direction, the system accelerates 
in that direction. This is equivalent to increasing the timestep in a steepest descent min¬ 
imization. Then, when the system overshoots the optimum, (V ■ 9””) < 0, the velocity is 
zeroed. 

This procedure is applied simultaneously to all images in the elastic band, /f and ff can 
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be derived from Cartesian components of the force nsing the following relations [56] : 


(F8) 

(F9) 


fi ^ fi cos Oi cos 4>i + ff cos Oi sin 4>i - U sin 
ft = - fi sin ipi + ff cos fi. 

When the VPO method is used to hnd an energy minimum, the force is simply the negative 
of the energy gradient with respect to spherical coordinates [SB]. The VPO optimization 
typically converges to the nearest local energy minimum. The global energy minimum can 
in principle be found by sampling enough initial guesses and repeating the procedure. 


2. Treatment of the poles 


During the force minimization, magnetic vectors may approach the vicinity of the poles 
where the azimuthal angle 0, becomes irrelevant. In this case the VPO formulated in terms 
of spherical coordinates will break down. Therefore, some other coordinate system needs 
to be chosen in the region of the poles. One option is to use the x and y components of 
magnetic moments. More specihcally, if < e or tt — < e, where e is some small positive 

number, then the coordinates mf and should be used instead of 9i and 0* in the VPO 
procedure: 


= sin^j cos0i, 
mf = sin^j sin0j. 


Furthermore, the velocity also needs to be expressed in terms of Cartesian coordinates [56] : 


nf = nf cos 6i cos 0* — vf sin 0*, 
Vf = cos 9i sin 0j + vf cos 0j. 


(FIO) 

(Fll) 


The evolution of a point in the e-neighborhood of the pole is then governed by the usual 
Newton’s equations of motion. The ff, ff and vf, vf components of the 2V-dimensional 
force and velocity vectors, 3^ and V, should be substituted by ff, ff and vf, vf in the 


velocity projection of Eq. (F7). When the system leaves the e-neighborhood of the pole. 


inverse transformation to spherical coordinates needs to be performed [56] . 

A typical trajectory of the representative point during the VPO calculations is shown in 
Fig. This trajectory has a cusp where the inner product of velocity and force becomes 
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negative and the velocity is zeroed. The minimization path changes direction at that point 
and subsequently the system advances slowly from at hrst but then accelerates as the force 
keeps pointing in a similar direction. The trajectory converges at the minimum which is 
located at the pole without a problem. 



FIG. 6. Minimization of the energy of a test system using the VPO algorithm as described in 
Appendix The trajectory has a cusp where the inner product of velocity and force becomes 
negative and the velocity is zeroed. The minimization path changes direction at that point and 
the system advances slowly from there at first but then accelerates as the force keeps pointing in 
a similar direction. The trajectory converges at the minimum which is located at the pole without 
a problem. 


Appendix G: Verification of a saddle point 

After a maximum along the MEP has been found, either by interpolation of the path 
between adjacent images or by using the CI-GNEB procedure, it is important to check 
whether this maximum indeed corresponds to a hrst order SP on the energy surface. This 
analysis is in any case needed to evaluate the transition rate using either HTST or Kramers 
methods. A hrst order SP is an extremum on the energy surface which is a maximum with 
respect to one and only one eigenmode while it is a minimum with respect to the other 
modes. In some cases, when a highly symmetrical initial path is used, the NEB can result in 
a path going through a higher order SP, i.e. a stationary points on the energy surface which 
is a maximum with respect to two or more modes. Energy surfaces can also be constructed 
where a ridge comes directly in continuation of an energy valley and the NEB method can 
then converge on a path lying partly along the ridge EH. For simple systems with only a few 
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degrees of freedom, as the ones considered in the test problems, the energy surface and the 
MEP can be visualized, so it is quite easy to check the results of the calculations. For complex 
magnetic systems with multiple magnetic vectors, a reliable strategy for verihcation of a SP 
is to calculate the 2N x 2N matrix of the second derivatives of the energy, the Hessian matrix, 
and calulate its eigenvalues. Any convenient set of coordinates can be used when computing 
the second derivatives, including spherical coordinates, stereographic coordinates, etc. If 
one and only one eigenvalue of the Hessian matrix is negative, then the maximum along the 
converged path is indeed a hrst order SP. If two eigenvalues are negative, then it is a second 
order saddle point, etc. If the NEB converges on a path going through a second order saddle 
point, then the highest energy image can be displaced in a direction along the negative 
mode that is orthogonal to the path tangent. Another option is to apply the minimum 
mode following method starting from the highest energy image 150] . Convergence to a 
path going through a higher order saddle point can also be the result of accidental symmetry 
in the initial path. It is in general good to add small random noise to the initial path to 
break such symmetries. 
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